Salmon Biphasic model 4

Author

Viktor Thunell

Published

September 4, 2025

Load libraries

Show code
# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "devtools","viridis","nls.multstart", "broom", "patchwork", "coda", "boot", "tidybayes","bayesplot", "nimbleHMC", "here")

# remotes::install_github("nimble-dev/nimble", ref = "devel", subdir = "packages/nimble")

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  }

invisible(lapply(pkgs, library, character.only = T))

options(ggplot2.continuous.colour = "viridis")
theme_set(theme_light()) # check why global theme option not working when rendering
# Set path
color_scheme_set("viridis")
home <- here::here()

Read data

Show code
sallaa <- readRDS(file = paste0(home,"/data/data-for-2-2/salmon-laa_2025-09-04.RData")) %>%
  rename(year = fi_year,
         site = sai_location,
         sea.age = sea_age_year,
         juv.age = juvenile_age_year,
         tot.age = tot_age_year,
         lat = fisa_y_4326,
         lon = fisa_x_4326) %>%
  #filter out non-aged individuals (~90000 individuals)
  filter(!(is.na(sea.age) & is.na(juv.age)))

Plot data

Show code
sallaa %>%
  drop_na(sex) %>%
  ggplot(aes(tot.age, length_mm, color = sex)) +
  geom_point() +
  facet_grid(age.type~sex) +
  scale_x_continuous(breaks = seq(0, 13, 1) ) +
  theme_light() 
drop_na: removed 125,641 rows (58%), 90,577 rows remaining

Show code
# If dropping all inds where both sea and sm age is NA
sallaa %>% 
  ggplot() +
  geom_density(aes(x = tot.age), fill = "deeppink", alpha = 0.5) +
  xlim(0, 13) +

sallaa %>% 
  ggplot() + 
  geom_density(aes(x = length_mm), fill = "deeppink", alpha = 0.5)

Show code
sallaa %>% 
  drop_na(sea.age) %>%
  ggplot() +
  geom_density(aes(x = length_mm, fill = factor(tot.age)), alpha = 0.5) +
  facet_wrap(~tot.age)
drop_na: removed 90,930 rows (42%), 125,288 rows remaining

Show code
# age at smoltification by spatial unit
sallaa %>% 
  filter(age.type == "both") %>%
  ggplot() +
  geom_density(aes(x = juv.age, color = spat.unit),  alpha = 0.5) +
  facet_wrap(~spat.unit) +
  xlim(-1, 5)
filter: removed 106,530 rows (49%), 109,688 rows remaining
Warning: Removed 2 rows containing non-finite outside the scale range
(`stat_density()`).

Show code
# juvenile size at age
sallaa %>%
  filter(age.type == "juve.only") %>%
  ggplot(aes(x = juv.age, y =length_mm, color = spat.unit)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = lm) +
  facet_wrap(~spat.unit)
filter: removed 125,288 rows (58%), 90,930 rows remaining
`geom_smooth()` using formula = 'y ~ x'

Show code
# size at age zero
sallaa %>%
  filter(!age.type == "sea.only",
         tot.age == 0) %>%
  ggplot() +
  geom_density(aes(x = length_mm), fill = "deeppink", alpha = 0.5) +
  labs(title = "Age 0")
filter: removed 165,288 rows (76%), 50,930 rows remaining

Biphasic 4

Spatiotemporal model for g and par.

Check - ages and sizes of smolts and pre smolts -

Show code
biph4.code <- nimbleCode({
  
  # likelihood
  for(i in 1:nobs){
    length_mm[i] ~ dnorm(l.mu[age.index[i], ind.id[i]], sd = sig.l)
    }
  
  # Estimate length at age for each individual (j in nn) and for each age of j (k in age.index=1 to age.index[j], i.e. age 0 to catch age)
  for(j in 1:nn){
    l.mu[1,j] <- lb.mu[spat.unit[j]]# old:+ g[spat.unit[j], hatch.year.f[j]] # age
    
    for(k in 2:age.index[j]+1){
    #Say that we for hatch.year=2010 want to predict length of an individual of age 3. If we then want to use environmental conditions in 2012 to predict length in 2013, year should be hatch.year + 2 = 2012 but here k = age + 1 = 4. Using k, we need hatch.year + k - 2. BUT, this makes for using hatch year for both k=1 (in l.mu[1,j]) and k=2. Furtherm. k-1 has different meaning in the step for smolt age.
    l.mu[k,j] <- l.mu[k-1,j] + step(smo.age[j] - (k-1))*exp(g[spat.unit[j], (hatch.year.f[j]+k-2)]) +
      (1-step(smo.age[j] - (k-1)))*((exp(par[spat.unit[j],(hatch.year.f[j]+k-2),1]) - l.mu[k-1,j])*(1 - exp(-exp(par[spat.unit[j],(hatch.year.f[j]+k-2),2]) #*(1 + (tot.age[j] %% 1))
      )))
    }
  }
  
  # estimate smo.age. 
  for(l in 1:nn){
    smo.age[l] ~ dgamma(shape = sh.sm[spat.unit[l]], scale = sc.sm)
    }
  
  # Priors  
  sig.l ~ dexp(1/500)
  #lb.mu ~ dnorm(17, sd = 5) # ~15 to 20 mm, https://doi.org/10.1111/j.1095-8649.2009.02497.x sd arbitary/guesstimate
  
  # LKJ prior on correlation matrix, see NIMBLE manual p45.
  Ustar[1:npars,1:npars] ~ dlkj_corr_cholesky(1.3, npars) # eta = 1.3
  U[1:npars,1:npars] <- uppertri_mult_diag(Ustar[1:npars, 1:npars], sig_par[1:npars])

  # priors for smolt age and juvenile growth rate
  
  sc.sm <- 0.5^2/2.5 # shape = 2.5^2/.5^2 = 25, #From mean=shape*scale, var=shape*scale^2 -> shape = mean^2/var & scale = var/mean. mean=2.5, sd=0.5.
  
  for(i in 1:nsu){
    sh.sm[i] ~ dnorm(25, sd = 5) # dgamma shape param, sd is arbitrary
    lb.mu[i] ~ dnorm(100, sd = 25) # FIX prior
    for(j in minyear:maxyear){ # Take into accuount hatch.year[j]+k-1
      par[i,j,1:npars] ~ dmnorm(mu_par[1:npars], cholesky = U[1:npars, 1:npars], prec_param = 0)
      g[i,j] ~ dnorm(mean = g_lmean, sd = g_lsd) # LOG SCALE and weak prior
    }
  }
  
  g_lsd <- sqrt(log(1 + (10^2) / (50^2))) # sqrt of variance (= sd) of the lognormal distr. 50 and 10 is arbitrary atm.
  g_lmean <- log(50) - 0.5 * g_lsd^2  # mean of the lognorm distr.
  
  mu_par[1] ~ dnorm(mean = l_lmean, sd = l_lsd)
  mu_par[2] ~ dnorm(mean = k_lmean, sd = k_lsd)
  sig_par[1] ~ dlnorm(0,1)
  sig_par[2] ~ dlnorm(0,1)

  # K & l values from fb  
  k_lsd <- sqrt(log(1 + (0.28^2) / (0.43^2))) 
  k_lmean <- log(0.43) - 0.5 * k_lsd^2  
  l_lsd <- sqrt(log(1 + (278^2) / (1378^2))) 
  l_lmean <- log(1378) - 0.5 * l_lsd^2
  })

# Function creating the Cholesky of the covar. matrix (p45 Nimble manual)
uppertri_mult_diag <- nimbleFunction(
  run = function(mat = double(2), vec = double(1)) {
    returnType(double(2))
    p <- length(vec)
    out <- matrix(nrow = p, ncol = p, init = FALSE)
    for(k in 1:p)
      out[ , k] <- mat[ , k] * vec[k]
    return(out)
   # turn off buildDerivs for the i index
}, buildDerivs = list(run = list(ignore = c('k'))))

data <- sallaa %>%
  filter(!age.type == "sea.only",
         !is.na(year),
         year > 2000
         ) %>% 
  slice_sample(n = 2000) %>%
  mutate(# create smolt age (!IMPROVE from data using fi_lfs_code)
         smo.age = if_else(age.type == "both", juv.age, NA),
         sea.age = if_else(is.na(sea.age), 0, sea.age),
         age.index = as.integer(tot.age + 1), # index for l.mu, + 1 to include age 0 individuals
         hatch.year = year-tot.age,
         hatch.year.f = as.numeric(factor(hatch.year, levels = sort(unique(hatch.year)), labels = seq_along(sort(unique(hatch.year))))),
         year.f = as.integer(factor(year)),
         spat.unit = as.integer(factor(spat.unit))) %>%
  select(smo.age, tot.age, length_mm, age.index, year, spat.unit, hatch.year, hatch.year.f) %>%
  mutate(ind.id = row_number()) # ind id
filter: removed 58,915 rows (27%), 157,303 rows remaining
slice_sample: removed 155,303 rows (99%), 2,000 rows remaining
mutate: changed 977 values (49%) of 'sea.age' (977 fewer NAs)
        converted 'spat.unit' from character to integer (0 new NA)
        new variable 'smo.age' (double) with 6 unique values and 49% NA
        new variable 'age.index' (integer) with 10 unique values and 0% NA
        new variable 'hatch.year' (double) with 31 unique values and 0% NA
        new variable 'hatch.year.f' (double) with 31 unique values and 0% NA
        new variable 'year.f' (integer) with 25 unique values and 0% NA
select: dropped 20 variables (sai_cou_code, site, origin, weight_g, sea.age, …)
mutate: new variable 'ind.id' (integer) with 2,000 unique values and 0% NA
Show code
npars <- 2
# get max year by, needed as hatch.year is integer and the last year is not (simpler solution?).  
maxyear <- max(data$hatch.year.f) + data %>% 
  filter(year == max(year)) %>% 
  filter(tot.age == max(tot.age)) %>% 
  distinct(tot.age) %>% pull(tot.age) 
filter: removed 1,997 rows (>99%), 3 rows remaining
filter: removed one row (33%), 2 rows remaining
distinct: removed one row (50%), one row remaining
Show code
# build model
biph4.model <- nimbleModel(biph4.code,
                           constants = list(npars = npars,
                                             nobs = nrow(data),
                                             nn = nrow(data), # change when >1 sample per ind.id
                                             nsu = length(unique(data$spat.unit)),
                                             age.index = data$age.index, 
                                             ind.id = data$ind.id,
                                             hatch.year.f = data$hatch.year.f,
                                             minyear = min(data$hatch.year.f),
                                             maxyear = maxyear,
                                             spat.unit = data$spat.unit),
                           data = data %>% select(-ind.id,-age.index,-spat.unit,-hatch.year,-tot.age, -year, -hatch.year.f))
select: dropped 7 variables (tot.age, age.index, year, spat.unit, hatch.year, …)
Defining model

Building model

Setting data and initial values

Running calculate on model
  [Note] Any error reports that follow may simply reflect missing values in model variables.

Checking model sizes and dimensions

  [Note] This model is not fully initialized. This is not an error.
         To see which variables are not initialized, use model$initializeInfo().
         For more information on model initialization, see help(modelInitialization).
Show code
# check values 
biph4.model$simulate("sig.l")
biph4.model$simulate("mu_par")
biph4.model$simulate("sig_par")
biph4.model$simulate("Ustar")
biph4.model$simulate("U")
biph4.model$simulate("sh.sm")
biph4.model$simulate("smo.age")
biph4.model$simulate("lb.mu")
biph4.model$simulate("g")
biph4.model$simulate("par")
biph4.model$simulate("l.mu")

biph4.model$sig.l
[1] 203.509
Show code
biph4.model$U
         [,1]      [,2]
[1,] 0.685916 0.2884687
[2,] 0.000000 0.6952038
Show code
biph4.model$Ustar
     [,1]      [,2]
[1,]    1 0.3832570
[2,]    0 0.9236417
Show code
biph4.model$mu_par
[1]  7.3382486 -0.3748789
Show code
biph4.model$sig_par
[1] 0.6859160 0.7526769
Show code
biph4.model$sh.sm
 [1] 22.80734 18.29895 24.84303 19.80169 22.10281 34.20938 28.30117 25.15342
 [9] 38.38293 29.50544 18.35988 29.36887
Show code
biph4.model$smo.age[1:50]
 [1] 1.000000 2.571189 1.000000 1.000000 2.780096 2.000000 1.000000 3.874723
 [9] 3.000000 2.000000 2.000000 2.440185 2.000000 2.301722 1.000000 2.702872
[17] 2.000000 2.884141 2.000000 1.000000 2.000000 2.000000 2.330448 3.000000
[25] 2.115186 1.794456 2.569748 3.000000 2.000000 2.000000 2.311045 4.100509
[33] 1.885151 2.394656 4.000000 1.000000 1.000000 1.000000 3.000000 2.757279
[41] 1.924156 3.000000 3.004745 3.000000 1.000000 2.501019 1.000000 3.546683
[49] 2.366453 1.000000
Show code
biph4.model$lb.mu
 [1] 104.75186 135.68843 110.01707 104.26566  90.84838 127.30670 101.52329
 [8]  98.82879 114.86248 144.76100 120.48958 102.77422
Show code
exp(biph4.model$g)
          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]
 [1,] 34.42720 51.72280 62.69168 58.11197 76.75703 60.17137 50.14673 41.98606
 [2,] 52.13390 69.49644 32.72035 43.07783 58.16661 39.93485 48.27357 59.55114
 [3,] 55.48917 41.66032 59.56387 45.86261 48.30688 45.91277 48.66094 48.44601
 [4,] 50.94406 44.67323 60.76316 46.21646 75.34816 47.26799 51.30011 57.79823
 [5,] 49.72428 46.87349 62.32715 38.18464 48.66287 67.17933 46.24349 56.46999
 [6,] 49.77200 77.36573 46.52932 46.45484 49.41016 51.75483 37.72351 58.45831
 [7,] 51.05327 49.97162 58.71447 47.88377 45.24087 63.90535 41.05441 59.83384
 [8,] 37.70689 55.28395 51.17743 54.90387 44.35532 44.89678 53.27905 55.31396
 [9,] 45.45734 52.86138 37.15070 64.11559 46.63762 45.56983 44.08551 46.80683
[10,] 43.08993 51.10170 78.11031 60.32844 54.33619 60.36043 73.04435 51.51648
[11,] 62.73698 48.83930 32.26505 46.95757 71.26317 38.77566 40.63142 44.20964
[12,] 36.93350 49.19554 44.92527 47.74869 44.55280 57.13482 40.24684 48.29059
          [,9]    [,10]    [,11]    [,12]    [,13]    [,14]    [,15]    [,16]
 [1,] 84.48650 54.86693 52.08447 48.76504 47.22146 44.78493 46.50656 52.67985
 [2,] 75.86255 99.59523 54.39426 61.83264 36.57759 38.97283 54.66969 57.64390
 [3,] 59.33554 42.34333 35.71948 54.55977 47.49285 38.68775 64.61671 52.13066
 [4,] 56.04016 57.67534 62.07442 48.39104 53.97144 45.85235 56.57665 38.23986
 [5,] 49.25813 44.19469 55.38866 63.20386 52.78787 49.50287 39.32551 42.84893
 [6,] 71.69852 61.82711 41.67357 38.11834 70.55941 44.91052 48.63787 41.16368
 [7,] 42.33127 51.33158 49.29930 43.86074 52.88007 56.76352 39.70058 45.62063
 [8,] 51.99418 59.86163 41.95414 46.68039 49.44480 42.60720 47.28234 38.32817
 [9,] 72.09801 53.28374 63.18165 51.42184 39.56409 44.93896 49.30640 45.48477
[10,] 48.43130 62.47986 44.81260 37.98589 41.84913 33.89100 39.36515 58.84307
[11,] 64.75171 55.87192 45.87807 47.99285 43.87252 44.36348 44.63569 47.87823
[12,] 54.48041 86.30090 35.84233 48.91697 54.10980 57.12564 41.83589 32.47545
         [,17]    [,18]    [,19]    [,20]    [,21]    [,22]    [,23]    [,24]
 [1,] 69.11628 56.54901 45.75586 34.47262 63.39292 34.41929 40.33396 54.51248
 [2,] 47.20973 45.67775 40.70438 52.05254 54.78839 56.60815 56.85560 64.13411
 [3,] 52.19172 50.96117 40.23503 52.12813 50.70336 52.39446 47.19573 43.77809
 [4,] 60.49847 65.38262 37.43927 44.01469 54.48135 52.74627 63.99180 48.53062
 [5,] 40.78562 44.87222 41.27387 69.51044 42.77725 75.69910 66.78379 36.81571
 [6,] 64.15869 44.10760 60.78312 56.13507 59.59881 61.74853 53.88143 41.13292
 [7,] 49.44624 48.91318 45.20188 59.94197 46.52275 33.41757 38.72768 48.15973
 [8,] 39.01122 39.33613 59.31604 46.45004 43.66240 47.32589 50.25836 55.05926
 [9,] 64.60066 58.16951 58.66229 78.92355 61.59924 61.57374 57.85964 34.58950
[10,] 40.43557 56.99610 53.73309 47.21610 49.27980 51.93589 59.72260 54.07465
[11,] 41.90699 60.02849 42.66656 37.81835 43.38051 44.40954 40.65779 44.13033
[12,] 72.25798 40.42952 41.79472 33.34564 53.68461 48.47118 63.89910 41.99203
         [,25]    [,26]    [,27]    [,28]    [,29]    [,30]    [,31]    [,32]
 [1,] 78.00881 46.18646 45.26350 53.44459 57.25141 46.18399 61.26333 50.19670
 [2,] 71.24540 47.20654 52.93293 53.36410 56.17169 52.36332 52.07289 48.72483
 [3,] 39.39310 84.28453 43.42337 44.13930 40.84783 33.74848 42.62128 43.94696
 [4,] 55.84752 62.44166 42.57960 51.59431 46.90590 53.35923 34.57439 55.19196
 [5,] 41.97635 46.86417 40.49910 60.61980 56.92274 32.12334 53.31458 46.28825
 [6,] 46.45075 48.86583 47.42269 68.12404 63.08481 45.95933 41.41478 61.40919
 [7,] 53.96756 56.77130 50.21567 58.52865 35.03986 57.45641 54.73099 56.03651
 [8,] 42.59742 50.40402 62.07492 77.05629 54.95880 64.27946 49.62921 34.46885
 [9,] 69.44832 57.07408 53.45280 40.24011 37.66502 45.71017 48.17530 38.47867
[10,] 40.23104 55.19771 39.28837 46.83396 45.95452 43.82402 42.45030 40.44739
[11,] 52.87033 52.89088 52.03459 47.33244 56.76787 74.58414 45.58892 43.71184
[12,] 57.63829 65.29669 46.51649 52.40497 55.80263 41.77943 44.91177 34.93486
         [,33]
 [1,] 39.66214
 [2,] 42.82060
 [3,] 47.44305
 [4,] 49.14023
 [5,] 45.04533
 [6,] 52.74020
 [7,] 53.48730
 [8,] 52.45705
 [9,] 42.65015
[10,] 47.42161
[11,] 53.15087
[12,] 42.91455
Show code
exp(biph4.model$par[1:11,1:20,1]) #linf
           [,1]      [,2]      [,3]      [,4]      [,5]     [,6]      [,7]
 [1,]  593.5666 3545.5926  685.3360  893.0542 1509.9457 4451.204 1602.3008
 [2,]  822.3401 1099.5602 2002.3381 1402.8689 1751.1519 1959.305 4351.6338
 [3,] 1265.7688 1127.8678 1128.9224 5490.3302 1313.7089 2276.209 1393.6186
 [4,] 1690.6379 2433.7213 1562.0013  429.8240 1888.9579 1780.323 1756.3152
 [5,] 1305.3765  574.0558 3106.7832  928.9854 1428.1928 1425.388  911.6412
 [6,] 2126.9635 1998.0317 1066.5719 1292.9898 2235.0225 1304.668 3305.5168
 [7,] 3128.9062 1178.2564 3459.5519 1685.5121 1317.1630 1210.917 1783.9054
 [8,]  596.1599 1662.5562  503.3115 1226.5548 1205.7703 3390.083 1961.5012
 [9,]  602.7349  820.7374 1824.2478 1766.8114  970.1561 2350.865  852.2550
[10,] 2033.3761 5100.4627 1161.6251  968.3964  801.9300 1011.133 1793.1289
[11,] 3377.9324  955.9493  570.9328  695.2532 1543.4400 3108.501  900.4184
           [,8]      [,9]     [,10]     [,11]     [,12]      [,13]     [,14]
 [1,] 1287.8216 2977.5120 1078.7364 3259.4207 1666.6069   417.9922 4121.4255
 [2,]  781.1206 3950.4005 1166.6360 1196.9912 2241.3951  1183.4963 2349.2522
 [3,] 1252.2748 1310.0631 1419.7015 2250.1357 3654.0039 12288.5433  372.2343
 [4,] 1118.0151 1861.4641  470.1928 1011.9702 1522.4650   611.7901 5478.1196
 [5,] 1534.8207 2062.8145 2886.7718 2726.1011  829.1002  2681.9950 1251.9641
 [6,] 2668.8411  826.5612 1756.9168  659.8585 1921.1699  1735.0257  928.4798
 [7,] 1120.2444 1732.3281 2102.7190 4340.2286 1824.0623  1292.0987 1290.4175
 [8,] 2367.4230 2648.5867 1302.6643 2563.5953 4839.7130  2687.7542 3880.1217
 [9,] 1195.6180 1025.8544  510.0036 3313.8378 1279.5249  2276.9070  673.1509
[10,] 3815.6908 1821.0199 1215.4660 1290.7135 3624.8442  2535.4602 4075.5236
[11,] 1440.6370 1772.4298 5262.0951 1409.6069 1676.5664  4179.9592  541.2000
          [,15]     [,16]     [,17]     [,18]     [,19]     [,20]
 [1,] 1047.9201 2489.0066 4338.6284 1434.7297  866.3850 1533.3274
 [2,] 1236.3297  438.3555 1897.5328 3193.8130  927.2877 2920.9706
 [3,] 1551.1936  970.3714 4855.4558  939.1186 2663.0389  955.4984
 [4,] 1356.0171 2602.5594 2949.3121 3625.7026 1288.3665 3029.4474
 [5,] 1647.4921 1984.2141 2201.7368 1425.2323 1545.2646  794.5688
 [6,] 2284.9896  973.1616 1356.3779 2117.6794 4828.4257 3314.8821
 [7,]  901.9795  510.5535  375.6534  463.7307 4230.5571 1246.1109
 [8,] 1335.6421  845.8185 1310.8643 1963.0006 1749.0603 1765.0289
 [9,]  633.3485  463.4456 1740.8399  943.3689 3339.4437 3862.9998
[10,] 1089.0547 3748.0046 1137.8435 1609.0560 3622.1472 2278.9221
[11,] 2711.3687 1371.3522 1327.8190 2492.6420 1411.0432  683.9628
Show code
exp(biph4.model$par[1:11,1:20,2]) #k
           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
 [1,] 0.3757857 0.6851028 0.6371013 0.4385733 0.6565294 1.1864357 0.3194203
 [2,] 0.6767621 0.5769388 0.1846687 1.1281815 0.3019901 0.5371295 0.9434316
 [3,] 0.4027349 0.7891651 0.6408985 0.5221832 0.3993972 0.7018097 0.1898211
 [4,] 3.0556308 0.7499321 0.4598650 0.2057216 1.6821092 0.7484576 0.7008861
 [5,] 0.6996827 0.6348633 0.5366471 0.4827800 1.4528715 0.3735826 0.3009558
 [6,] 0.6832037 1.1421346 0.6879146 0.5429381 1.0459664 0.4632092 1.3335233
 [7,] 1.4367357 0.5352937 0.4354562 0.4198319 1.3763749 0.6329482 1.2931431
 [8,] 1.0548506 0.3390614 0.3081547 0.4092549 1.3336851 1.6922207 0.8876120
 [9,] 0.4656071 0.5153398 1.1538191 1.8472528 0.5793302 0.1722491 0.9320798
[10,] 1.0324939 0.5407047 0.5034960 0.5801327 0.2670956 0.4807999 1.3222879
[11,] 0.4036233 0.2682479 0.5358320 0.3438721 1.8083353 0.7908408 0.3105159
           [,8]      [,9]     [,10]     [,11]     [,12]     [,13]     [,14]
 [1,] 0.6930383 1.8522382 0.5657771 0.6264887 0.3905470 0.6207832 2.6947077
 [2,] 0.3219206 0.7413521 0.3515261 1.9397310 1.2212404 0.8221683 1.6406807
 [3,] 0.5414938 0.4505846 2.6754597 0.3210424 1.3087026 4.0714782 0.4439580
 [4,] 0.4476096 0.8007253 0.2399298 0.3188475 2.5158172 0.2146457 0.6590223
 [5,] 2.7075889 2.3785053 1.1213890 0.5966439 0.3246170 0.8297732 0.2657363
 [6,] 0.3418760 0.4904718 0.8603294 1.2589028 1.4885692 1.1657748 0.6045861
 [7,] 0.9411432 0.4365831 0.6748903 2.2057390 0.2292631 0.5337172 0.1634621
 [8,] 2.6084498 1.7111780 1.5190758 0.6866697 1.0542213 0.7796873 0.3873415
 [9,] 0.4994036 1.3256530 0.2095832 1.1704239 1.1851398 0.8786183 2.8939298
[10,] 0.1874818 1.5053492 1.2630751 0.9256293 0.9466119 0.2435952 1.6924856
[11,] 1.3072500 0.1102466 1.8639102 0.5035988 0.7777301 0.3453607 0.5164525
          [,15]     [,16]     [,17]     [,18]     [,19]     [,20]
 [1,] 0.4624555 0.7544328 1.7053628 1.3125130 0.3937441 0.6626281
 [2,] 0.4386369 1.4178890 0.7787941 2.0701606 0.3186812 0.4323573
 [3,] 0.7011044 0.5194709 0.5141824 0.4476104 0.5787509 0.7191883
 [4,] 0.4237093 1.3482122 0.9142583 1.6338722 1.5677129 1.1505732
 [5,] 0.7284173 0.7027967 0.7863591 0.6870926 0.8742686 0.6535705
 [6,] 0.7933637 0.5311352 0.2478311 0.6180635 1.1510170 0.9151332
 [7,] 0.4330051 0.3924334 0.3306910 0.6974283 4.5795571 0.7826731
 [8,] 0.2697981 0.1865165 0.8765391 1.2697006 0.4123909 0.7841659
 [9,] 0.7515729 0.4620418 0.6210606 0.2486285 1.6606960 1.6695201
[10,] 0.2869051 0.6678507 0.4075714 2.6994471 1.2031784 0.6988340
[11,] 0.3878411 1.2071617 0.4457232 0.7194973 0.9718995 0.5891420
Show code
# NA l.mus
sum(is.na(biph4.model$l.mu))/length(biph4.model$l.mu) # percent negatives
[1] 0.6297273
Show code
# NAs are after
biph4.model$l.mu[,1:10] # first 10 individuals
           [,1]     [,2]       [,3]       [,4]      [,5]     [,6]      [,7]
 [1,]  104.7519 144.7610   98.82879   98.82879  98.82879 114.8625  144.7610
 [2,]  153.5169 191.9771  137.15697  137.84001 149.23281 154.4266  196.6969
 [3,]  275.8306       NA  822.34277 1450.28433        NA 199.3655  996.7027
 [4,] 3861.6090       NA 1642.57191 1551.25106        NA 428.6716 5852.4332
 [5,]        NA       NA 1678.55802 1667.43926        NA 441.5382        NA
 [6,]        NA       NA 1725.55493 1514.24446        NA       NA        NA
 [7,]        NA       NA         NA 1396.26817        NA       NA        NA
 [8,]        NA       NA         NA         NA        NA       NA        NA
 [9,]        NA       NA         NA         NA        NA       NA        NA
[10,]        NA       NA         NA         NA        NA       NA        NA
[11,]        NA       NA         NA         NA        NA       NA        NA
          [,8]      [,9]      [,10]
 [1,] 144.7610  101.5233   98.82879
 [2,] 196.6969  165.4286  140.78293
 [3,]       NA  206.4830  187.46331
 [4,]       NA  266.3169 1541.24737
 [5,]       NA  784.9319 2292.35538
 [6,]       NA 1431.6856         NA
 [7,]       NA        NA         NA
 [8,]       NA        NA         NA
 [9,]       NA        NA         NA
[10,]       NA        NA         NA
[11,]       NA        NA         NA

Sample g, par, lb.mu, sig.l

Show code
# configure..
biph4.confmcmc <- configureMCMC(biph4.model, monitors = c("g","par","lb.mu","sig.l"))
===== Monitors =====
thin = 1: g, lb.mu, par, sig.l
===== Samplers =====
RW_lkj_corr_cholesky sampler (1)
  - Ustar[1:2, 1:2] 
RW_block sampler (266)
  - par[]  (266 multivariate elements)
RW sampler (656)
  - sig.l
  - sh.sm[]  (12 elements)
  - sig_par[]  (2 elements)
  - smo.age[]  (373 elements)
  - g[]  (266 elements)
  - mu_par[]  (2 elements)
posterior_predictive sampler (864)
  - smo.age[]  (604 elements)
  - g[]  (130 elements)
  - par[]  (130 multivariate elements)
conjugate sampler (12)
  - lb.mu[]  (12 elements)
Show code
# and build MCMC
biph4.mcmc <- buildMCMC(biph4.confmcmc)

# compile model
biph4.c <- compileNimble(biph4.model)
Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Show code
# compile mcmc (and specify the project model)
biph4.mcmcc <- compileNimble(biph4.mcmc)
Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Show code
# MCMC Samples
biph4.samp <- runMCMC(biph4.mcmcc, niter = 20000, nburnin = 15000, nchains = 2, samplesAsCodaMCMC = TRUE)
running chain 1...
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
running chain 2...
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

Sample l.mu

Show code
# configure..
biph4.confmcmc <- configureMCMC(biph4.model, monitors = c("l.mu"))
===== Monitors =====
thin = 1: l.mu
===== Samplers =====
RW_lkj_corr_cholesky sampler (1)
  - Ustar[1:2, 1:2] 
RW_block sampler (266)
  - par[]  (266 multivariate elements)
RW sampler (656)
  - sig.l
  - sh.sm[]  (12 elements)
  - sig_par[]  (2 elements)
  - smo.age[]  (373 elements)
  - g[]  (266 elements)
  - mu_par[]  (2 elements)
posterior_predictive sampler (864)
  - smo.age[]  (604 elements)
  - g[]  (130 elements)
  - par[]  (130 multivariate elements)
conjugate sampler (12)
  - lb.mu[]  (12 elements)
Show code
# and build MCMC
biph4.mcmc.lmu <- buildMCMC(biph4.confmcmc)

# compile model
biph4.c <- compileNimble(biph4.model)
Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Show code
# compile mcmc (and specify the project model)
biph4.mcmcc.lmu <- compileNimble(biph4.mcmc.lmu)
Compiling
  [Note] This may take a minute.
  [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Show code
# MCMC Samples
biph4.samp.lmu <- runMCMC(biph4.mcmcc.lmu, niter = 15000, nburnin = 10000, nchains = 2, samplesAsCodaMCMC = TRUE)
running chain 1...
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|
running chain 2...
|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|

Check g, par, lb.mu, sig.l

Show code
samp <- biph4.samp
sum <- summary(samp)

# summary of values of sampled parameters
vars_subset <- grep("^g\\[", rownames(sum$statistics), value = TRUE)
summary(exp(sum$statistics[which(rownames(sum$statistics)%in%vars_subset),][,1]))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  41.63   48.05   48.97   48.46   49.14   53.93 
Show code
vars_subset <- grep("^par\\[", rownames(sum$statistics), value = TRUE)
vars_subset1 <- grep("1]", vars_subset, value = TRUE)
summary(exp(sum$statistics[which(rownames(sum$statistics)%in%vars_subset1),][,1]))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  722.9   974.1   980.9   984.3   984.7  1429.1 
Show code
vars_subset2 <- grep("2]", vars_subset, value = TRUE)
summary(exp(sum$statistics[which(rownames(sum$statistics)%in%vars_subset2),][,1]))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3455  0.8524  0.8611  0.8812  0.8771  1.6282 
Show code
vars_subset <- grep("^lb.mu\\[", rownames(sum$statistics), value = TRUE)
summary(sum$statistics[which(rownames(sum$statistics)%in%vars_subset),][,1])
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  19.57   63.67   72.19   72.16   84.12  121.13 
Show code
vars_subset <- grep("sig.l", rownames(sum$statistics), value = TRUE)
sum$statistics[which(rownames(sum$statistics)%in%vars_subset),]
          Mean             SD       Naive SE Time-series SE 
  38.651434433    0.662647114    0.006626471    0.018106568 
Show code
# some trace plots
rn <- sample(1:200, 20)
vars_subset <- grep("^g\\[", varnames(samp), value = TRUE)[rn]
traceplot(samp[, vars_subset, drop = FALSE])

Show code
vars_subset <- grep("^par\\[", varnames(samp), value = TRUE)[rn]
traceplot(samp[, vars_subset, drop = FALSE])

Show code
vars_subset <- grep("^lb.mu\\[", varnames(samp), value = TRUE)
traceplot(samp[, vars_subset, drop = FALSE])

Show code
vars_subset <- grep("sig.l", varnames(samp), value = TRUE)
traceplot(samp[, vars_subset, drop = FALSE])

Plot g, par…

Show code
# g[spat.unit, hatch.year]
samp %>%
  gather_draws(g[su, hy], sep = ",") %>%
  filter(hy == c(1:5)) %>%
  ggplot() +
  geom_density(aes(x = .value, color = factor(hy))) + # age parameter
  facet_wrap(~su, scales = "free") +
  theme_light()
filter (grouped): removed 3,840,000 rows (97%), 120,000 rows remaining (removed
336 groups, 60 groups remaining)

Show code
# lb.mu[spat.unit]
samp %>%
  gather_draws(lb.mu[su], sep = ",") %>%
  #filter(hy == c(1:5)) %>%
  ggplot() +
  geom_density(aes(x = .value, color = factor(su))) + # age parameter
  #facet_wrap(~su, scales = "free") +
  theme_light()

Show code
# growth curves not looking ok. Should plot the individual curves and summarise the mean by su curve
samp %>%
  spread_draws(par[su,hy,p], g[su,hy], lb.mu[su], sep = ",") %>%
  mutate(g = exp(g),
         par = exp(par)) %>%
  group_by(.draw,su,p) %>%
  summarise(par.m = mean(par), g.m = mean(g), lb.mu = mean(lb.mu), .groups = "keep") %>%
  group_by(su,p) %>%
  median_qi() %>%
  mutate(p = case_when(p == 1 ~ "Linf",
                       p == 2 ~ "K",
                       .default = NA)) %>%
  pivot_wider(id_cols = su, names_from = p, values_from = c(par.m, g.m, lb.mu)) %>%
  expand_grid(age = seq(0, 12, 0.1)) %>%
  mutate(length.a = par.m_Linf*(1-exp(-par.m_K*age)),
         length.j = lb.mu_Linf + g.m_Linf) %>%
  ggplot() +
  geom_point(data = sallaa %>% filter(!age.type == "sea.only"), 
             aes(tot.age, length_mm), alpha = 0.1, color = "lightblue") +
  geom_line(aes(age, length.a, color = factor(su))) +
  geom_line(aes(age, length.j, color = factor(su))) +
  theme_light() 
mutate (grouped): changed 7,920,000 values (100%) of 'par' (0 new NAs)
                  changed 7,920,000 values (100%) of 'g' (0 new NAs)
group_by: 3 grouping variables (.draw, su, p)
summarise: now 240,000 rows and 6 columns, 3 group variables remaining (.draw, su, p)
group_by: 2 grouping variables (su, p)
mutate: converted 'p' from integer to character (0 new NA)
pivot_wider: reorganized (p, par.m, par.m.lower, par.m.upper, g.m, …) into (par.m_Linf, par.m_K, g.m_Linf, g.m_K, lb.mu_Linf, …) [was 24x14, now 12x7]
mutate: new variable 'length.a' (double) with 1,441 unique values and 0% NA
        new variable 'length.j' (double) with 12 unique values and 0% NA
filter: removed 15,600 rows (7%), 200,618 rows remaining

Check l.mu

Show code
#samp.lmu <- biph4.samp.lmu
# remove NA nodes (the older-than-catch-age nodes)
cm <- lapply(biph4.samp.lmu, as.matrix)
nonNA_nodes <- cm %>% 
  map(~ colnames(.x)[colSums(is.na(.x)) == 0]) %>%
  reduce(intersect) # applies function intersect over elements of the lists returned by map(), i.e. removing the colnames from both lists.
cm_cleaned <- lapply(cm, function(x) x[, nonNA_nodes, drop = FALSE])
samp.cl <- mcmc.list(lapply(cm_cleaned, mcmc))  

#sum <- summary(samp.cl) # not working when l.mu == NA

# some traces of l.mu
vars_subset <- grep("^l.mu\\[", varnames(samp.cl), value = TRUE)[rn]
traceplot(samp.cl[, vars_subset, drop = FALSE])

Plot predicted vs observed l.mu

Show code
# plot observed against predicted
obs <- data %>%
  select(length_mm,tot.age,ind.id) %>%
  mutate(tot.age = round(tot.age),
         observed = length_mm)
select: dropped 6 variables (smo.age, age.index, year, spat.unit, hatch.year, …)
mutate: new variable 'observed' (double) with 392 unique values and 0% NA
Show code
predi <- samp.cl %>%
  gather_draws(l.mu[age,ind.id], sep = ",") %>%
  median_qi() %>%
  rename(predicted = .value) %>%
  mutate(tot.age = age-1) %>%
  filter(tot.age %in% obs$tot.age & ind.id %in% obs$ind.id)
rename: renamed one variable (predicted)
mutate: new variable 'tot.age' (double) with 11 unique values and 0% NA
filter: removed 2 rows (<1%), 8,144 rows remaining
Show code
predi %>%
  left_join(obs, by = join_by(tot.age, ind.id)) %>%
  ggplot(aes(x = observed, y = predicted)) +
  geom_point(size = 2, alpha = 0.7) +
  #geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  theme_minimal()
left_join: added 2 columns (length_mm, observed)
           > rows only in x    6,144
           > rows only in obs (    0)
           > matched rows      2,000
           >                  =======
           > rows total        8,144
Warning: Removed 6144 rows containing missing values or values outside the scale range
(`geom_point()`).

Show code
# smolt ages
# data %>%
#   filter(!is.na(smo.age)) %>%
#   select(smo.age) %>%
#   mutate(t= "obs") %>%
#   bind_rows( samp.cl %>% gather_draws(smo.age[i]) %>% rename(smo.age = .value) %>% mutate(t="pred") ) %>%
#   ggplot() +
#   geom_density(aes(smo.age, fill = t))